%==========================================================================
%=  This file is part of the Sound Restoration Project
%=  (c) Copyright Industrial Mathematics Institute
%=                University of South Carolina, Department of Mathematics
%=  ALL RIGHTS RESERVED
%=
%=  Author: Borislav Karaivanov
%==========================================================================

%==========================================================================
% List of the files on which this procedure depends:
%
% fitReduceFitSignal.m
% Curve Fitting toolbox
%
%==========================================================================

%==========================================================================
% The function "findBestFitToSignal" takes an array of signals of the same
% length and calibrates each of them with an individual calibration mask
% based on the average of a given number of signals closest to the one
% being calibrated. For the signals far from the beginning and the end of
% the array of signals the signals entering the mask are positioned
% symmetrically on both sides of the calibrated signal while for the
% signals near the beginning and the end the masks still use the same
% number of closest signals but those are asymmetrically positioned as one
% side of the calibrated signal there will be shortage of close signals.
% INPUT: "signalArr" is a 2D array whose columns hold the individual
% signals to be calibrated.
% "numSmoothings" is an optional non-negative integer specifying the number
% of smoothing passes to be applied to the given signal.
% "movingAveSpan" is an optional, positive, odd integer specifying the span
% of the stencil to be used for moving average smoothing of the given
% signal.
% "degOfPolyFitArr" is an optional array of non-negative integers
% specifying the degrees of the polynomials to be fitted to the given
% signal. If this parameter is specified (and non-empty), then the values
% of the previous two parameters are ignored.
% "maxAllowedCut" is an optional, non-negative integer specifying the
% maximal allowed cumulative number of samples that can be removed at the
% beginning or at the end of the given signal in order to produce a reduced
% signal with a better polynomial fit.
% OUTPUT: "fitArr" returns a vector of length not greater than that of the
% original signal holding the best fit to the reduced version of the given
% signal.
% "leftCut" returns the number of samples removed at the beginning of the
% given signal.
% "rightCut" returns the number of samples removed at the end of the given
% signal.
% "degOfPolyFit" returns the degree of polynomial for which the best fit is
% attained, if polynomial fits were computed, or an empty array ([]) if the
% signal was smoothed.
%==========================================================================
function [fitArr, leftCut, rightCut, degOfPolyFit] = findBestFitToSignal(signalArr, ...
    numSmoothings, movingAveSpan, degOfPolyFitArr, maxAllowedCut)

if (nargin < 2)
    numSmoothings = 1;
end
if (nargin < 3)
    movingAveSpan = 15;
end
if (nargin < 4)
    degOfPolyFitArr = (1:4);
end
if (nargin < 5)
    maxAllowedCut = 100;
end


% Take different courses depending on whether the given signal is to be
% smoothed or fitted.
if ( (nargin < 4)||(isempty(degOfPolyFitArr) == true) )
    
    % Smooth the given signal and its reduced versions by applying the
    % desired number of moving average passes.
    fitFunctHandle = @(x) repeatedMovingAverage(x, numSmoothings, movingAveSpan);
    [fitArr, leftCut, rightCut] = ...
        fitReduceFitSignal(signalArr, fitFunctHandle, maxAllowedCut);
    
    degOfPolyFit = [];
    
else
    
    % Allocate memory for the results from the different polynomial fits.
    fitCellArr = cell(size(degOfPolyFitArr));
    maxDiffArr = zeros(size(degOfPolyFitArr));
    leftCutArr = zeros(size(degOfPolyFitArr));
    rightCutArr = zeros(size(degOfPolyFitArr));
    
    for k = 1:numel(degOfPolyFitArr)
        
        % Set a handle to the function implementing the polynomial fit of
        % the current degree.
        fitFunctHandle = @(x) polynomialFit(x, degOfPolyFitArr(k));
        
        % Compute the desired polynomial fit to the reduced signal.
        [fitCellArr{k}, leftCutArr(k), rightCutArr(k), maxDiffArr(k)] = ...
            fitReduceFitSignal(signalArr, fitFunctHandle, maxAllowedCut);
        
    end
    
    % Find the best of the computed polynomial fits.
    [~, minMaxDiffInd] = min(maxDiffArr);
    fitArr = fitCellArr{minMaxDiffInd};
    leftCut = leftCutArr(minMaxDiffInd);
    rightCut = rightCutArr(minMaxDiffInd);
    degOfPolyFit = degOfPolyFitArr(minMaxDiffInd);
end

return;
% end of the function "findBestFitToSignal"



%==========================================================================
% The function "polynomialFit" smoothes a given signal by applying
% a specified number of moving average passes.
% INPUT: "signalArr" is a vector holding the signal to be smoothed.
% "numSmoothings" is a non-negative integer specifying the number of
% smoothing passes to be applied to the given signal.
% "movingAveSpan" is a positive, odd integer specifying the span of the 
% stencil to be used for moving average smoothing of the given signal.
% OUTPUT: "signalArr" returns a vector holding the smoothed signal.
%==========================================================================
function signalArr = polynomialFit(signalArr, degOfPolyFit)

% Fit the given signal with a polynomial of the desired degree.
xArr = (1:numel(signalArr)).';
curveFunct = fit(xArr, signalArr, ['poly' num2str(degOfPolyFit)]);

% Sample the found polynomial.
signalArr = curveFunct(xArr);

return;
% end of the function "polynomialFit"



%==========================================================================
% The function "repeatedMovingAverage" smoothes a given signal by applying
% a specified number of moving average passes.
% INPUT: "signalArr" is a vector holding the signal to be smoothed.
% "numSmoothings" is a non-negative integer specifying the number of
% smoothing passes to be applied to the given signal.
% "movingAveSpan" is a positive, odd integer specifying the span of the 
% stencil to be used for moving average smoothing of the given signal.
% OUTPUT: "signalArr" returns a vector holding the smoothed signal.
%==========================================================================
function signalArr = repeatedMovingAverage(signalArr, numSmoothings, movingAveSpan)

% Smooth signal by applying the desired number of moving average passes.
for k = 1:numSmoothings
    signalArr = smooth(signalArr, movingAveSpan, 'moving');
end

return;
% end of the function "repeatedMovingAverage"
